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ABSTRACT 

This  paper  first  describes  a modified  version  of  the 
parametric  principal  pivoting  algorithm  [4]  for  the  class  of 
parametric  linear  complementarity  problems  with  P-matriccs. 

The  modified  version  of  the  parametric  principal  pivoting  al- 
gorithm is  then  used  to  develop  a new  and  efficient  algorithm 
to  solve  the  class  of  portfolio  analysis  problems  with  positive 
definite  (or  equivalently,  nonsingular)  covariance  matrices. 
Extension  of  the  new  algorithm  to  handle  explicit  upper-bounds 
is  also  established.  The  new  algorithm  and  its  extension  are 
then  specialized  to  the  "index  models"  introduced  by  Sharpe 
l28l  / I29]  • In  these  specializations,  the  algorithms  are  par- 

ticularly effective,  achieving  dramatic  savings  in  both  storage 
and  computations. 
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SOME  NEW  AND  EFFICIENT  ALGORITHMS 
FOR  PORTFOLIO  ANALYSIS 


Jong-Shi  Pang 


1.  INTRODUCTION 

In  1/952,  Markowitz  [19]  suggested  formulating  the  general  portfolio 
analysis  problem  as  a parametric  (convex)  quadratic  programming  problem. 
The  quadratic  term  in  the  objective  function  corresponds  to  the  predicted 
variance  of  portfolio  return  which  is  used  as  a measure  of  risk;  the  lin- 
ear term  in  the  objective  function  corresponds  to  the  expected  value  of 
portfolio  return;  the  parameter  of  the  problem  which  is  assumed  to  be 
nonnegative  and  appears  in  the  linear  term  of  the  objective  is  a measure 
of  the  relative  importance  of  reduced  risk  and  increased  return;  finally, 
the  constraints  consist  of  the  nonnegativity  of  the  proportions  of  an  in- 
vestor's fund  to  be  invested  in  the  securities  and  the  sum  of  the  propor- 
tions being  equal  to  unity.  Specifically,  the  general  portfolio  analysis 
problem  is 


U.l) 


minimize  tt  V - 0E 

2 p p 

n 

with  V - Var(y  x.R.) 

P i=l  " ^ 


n 

Ep  = r.xpil 


subject  to 
and 


i=l 
X > 0 


T 1 

c X = 1 


where  n is  the  number  of  securities;  0 is  the  nonnegative  parameter; 

T n 

X = (Xj^,...,x^)  e R is  the  vector  of  proportions  invested  in  the  secur- 

T 


ities;  R = (R^ , . . . ,R^) 


is  the  vector  of  returns  on  the  securities 
which  are  assumed  to  be  stochastic;  and  c is  the  n-vector  of  I's  . 
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Letting  V be  the  covariance  matrix  of  the  security  returns,  we  can  re- 
write problem  (1.1)  in  the  following  equivalent  format. 


(1.2) 


. . IT  T 

minimize  j x Vx  + 0r  x 


subject  to  X ^ 0 
T 

and  e X = 1 


where  r = (r^)  c R is  the  vector  with  r^  = -Exp(R^)  for  every 

i = 1, . . . ,n  . 

The  matrix  V,  being  a covariance  matrix,  is  symmetric  and  positive 
semi -dcf ini te . Therefore,  for  every  fixed  9,  problem  (1.2)  is  a convex 
quadratic  program.  Methods  for  solving  (1.2)  parametrically  (i.e.  obtain- 
ing a solution  for  every  value  of  6)  include,  among  others,  Marl^owitz' 
critical  line  method  (20)  and  Wolfe's  simplex-li):e  parametric  algorithm 
|3S).  Mar)\owitz'  method  starts  by  finding  a portfolio  with  maximum  ex- 
pected return  (i.e.  a solution  corresponding  to  0 = °=) , then  decreases  the 
Value  of  0 successively  and  proceeds  in  parametric  manner.  The  algorithm 

when  it  reaches  the  minimum  risk  portfolio  (i.e.  when  6 = 0).  Wolfe's 
parametric  algorithm  operates  in  precisely  the  reverse  order.  It  starts 
by  using  Uie  unparametarized  version  of  the  algorithm  to  obtain  an  initial 
solution  corresponding  to  0 = 0 and  then  increases  the  value  of  6 until 
it  reaches  Both  methods  (and  a few  others)  require  the  storage  of  the 

entire  matrix  V and  an  extensive  amount  of  computational  effort. 

Recognizing  the  informational  and  computational  complexity  of  the 
general  model,  Sharpe  (28)  introduced  some  simplified  models  for  portfolio 
analysis.  These  are  the  "index  models"  which  assume  that  the  security  re- 
turns arc  expressed  in  terms  of  some  market  indices.  Specifically,  a 
typical  index  model  assumes  that 


(1.3) 


R.  =ci.  +B,.I,  + ...+P.I  + c. 

a i 111  mim  i 


i = 1, . . . ,n 
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where  the  a's  and  g's  are  constant  scalars,  the  I's  are  appropriate 
market  indices  which  are  random  variables  with  mean  zero,  and  the  e's  are 
random  variables  with  mean  zero.  There  are  three  assumptions  to  the  model: 


(1.4a) 

Cov(I^, 

^3’ 

= 0 

for 

i + j 

(1.4b) 

Cov(c^, 

"j' 

= 0 

for 

i + j 

(1.4c) 

Cov ( c ^ , 

^3’ 

= 0 

for 

all  i,j 

Under  these  assumptions,  as  shown  in  [30 ] , the  covariance  matrix  V is 
given  by 


111  m 

(1.5)  V = 5;  + 5;  Var(I  .)  p 

j=l  333 

where  ^ is  a nonnegative  diagonal  matrix  whose  i-th  diagonal  entry  is 
equal  to  Var(c^),  and  for  j = l,...,m  . If  m = 1, 

the  model  is  called  a single  index  model , and  in  general,  it  is  called  an 
m-index  model.  P’or  practical  purposes,  it  is  reasonable  to  assume  that  m 
is  smaller  than  n . 

One  of  the  simplifications  resulting  from  the  introduction  of  the 
index  models  was  the  transformation  of  the  quadratic  term  x"^Vx  in  the 
objective  function  of  (1.2)  into  one  with  a diagonal  matrix  (i.e.  into  a 
sum  of  square  terms).  By  (1.5),  problem  (1.2)  can  be  cast  in  the  form: 

(1.6)  minimize  | (J  J)  {p  + GC^)"^  (^) 


subject  to 


where  y = <yi'  • • • 'Vni^ 
the  equalities 


X > 0 


e'^x  = 1 


yj  = Var(Ij)  PjX  j = l,...,m  . 

By  treating  the  y's  as  separate  variables  and 


y^  = Var(Ij)  p^x 


j = l,...,m 
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as  additional  constraints,  special  purpose  quadratic  programming  codes  have 
been  written  [12],  [34l  to  take  full  advantage  of  the  fact  that  the  objec- 
tive now  has  a diagonal  matrix. 

Sharpe  [31]  and  Stone  [32]  have  proposed  another  approach  to  simplify 
the  general  model.  This  involves  first  diagonalizing  the  covariance  matrix 
V by  some  numerical  procedure,  e.g.  the  one  described  in  Dantzig  [7, 
p.  491-492] , then  approximating  the  resulting  diagonalized  objective  func- 
tion by  a piecewise  linear  function,  and  finally,  applying  linear  program- 
ming to  solve  the  approximated  problem.  With  the  index  models,  the  first 
step  of  diagonalizing  the  covariance  matrix  V is  not  required  because  of 
the  formulation  (].6).  Of  course,  the  solution  obtained  by  this  linear 
programming  approach  gives  only  an  approximated  solution  to  the  general 
problem  (1.2) . 

Due  to  some  legal  or  personal  constraints,  it  is  common  in  practice 
to  impose  explicit  upper  bounds  on  the  proportions  invested  in  the  secur- 
ities. An  example  of  this  was  described  in  [29]  concerning  the  mutual 
fund  portfolio  analysis.  To  formulate  the  general  portfolio  analysis 
problem  with  explicit  upper  bounds,  let  a be  a positive  vector,  then  the 
problem  is 

It  T 

(1.7)  minimize  j ^ Vx  + Or  x 

subject  to  a ^ X ^ 0 

T 

and  e X = 1 

Assuming  a single  index  model  and  that  the  bounds  are  rather  restric- 
tive, Sharpe  [29]  proposed  a linear  programming  approximation  to  problem 

(1.7) .  Although  the  solution  procedure  then  becomes  a very  simple  one, 

(as  a matter  of  fact,  the  approximated  problem  is  almost  trivial)  it  pro- 
duces only  an  approximate  solution  to  the  original  problem. 

Recently,  Jucker  and  de  Faro  [13],  assuming  that  the  covariance 
matrix  V is  diagonal,  described  a simple  algorithm  for  problem  (1.2)  for 
fixed  value  of  the  parameter  6 ind  extended  the  algorithm  to  the  problem 


(1.7)  with  upper  bounds.  They  also  provided  a graphical  interpretation  to 
their  algorithm.  As  a matter  of  fact,  the  problem  they  considered,  namely, 
(1.2)  with  diagonal  V and  fixed  6 , is  a special  case  of  the  single  index 
model  corresponding  to  Var(Ij^)  = 0 . More  recently,  Alexander  [1]  has  re- 
ported some  computational  experience  with  the  use  of  Lem)<.e*s  algorithm  [17] 
to  solve  problem  (1.2)  in  a somewhat  different  formulation.  The  largest 
problem  solved  in  the  paper  contains  90  securities  and  was  solved  in  an 
average  of  190.3  seconds  (CPU  time).  Finally,  Rosenberg  and  Rudd  [24]  has 
reported  some  rather  encouraging  computational  results  with  the  use  of  a 
new  and  general-purpose  quadratic  programming  algorithm  in  portfolio  anal- 
ysis. The  algorithm  is  due  to  Von  Hohenbalken  [33]  and  is  not  parametric. 

It  uses  techniques  in  linear  programming. 

Despite  tlie  fact  that  there  are  many  quadratic  programming  algorithms 
that  can  be  used  to  solve  the  problems  (1.2)  and  (1.7),  few  of  them  are 
successful  in  solving  problems  with  a large  number  of  securities  (say, 
when  n is  in  the  order  of  hundreds  or  thousands).  A major  limitation  in 
the  application  of  tliese  algorithms  to  large-scale  problems  is  the  exces- 
sive computational  cost.  This  is  very  often  caused  by  the  lack  of  effici- 
ency in  the  implementation  of  the  algorithms. 

In  this  paper,  we  propose  a fresh  approach  and  develop  some  new  al- 
gorithms to  solve  the  classes  of  portfolio  problems  (1.2)  and  (1.7)  with 
positive  definite  covariance  matrices  V . We  will  also  discuss  efficient 
implementation  of  these  algoritlims,  with  special  emphasis  placed  on  tlieir 
implementation  for  large-scale  problems.  It  sliould  be  pointed  out  that 
since  V is  a covariance  matrix,  it  is  positive  definite  if  and  only  if 
it  is  nonsingular.  The  essential  tool  used  to  develop  the  new  algorithms 
is  a modified  version  of  the  parametric  principal  pivoting  algorithm  (4) 
for  the  class  of  parametric  linear  complementarity  problems  with  P-matrices. 
We  sliall  explain  the  terminology  and  describe  the  modified  algorithm  in  the 
next  section. 

Tl)e  new  algorithms  that  we  shall  propose  for  the  general  portfolio 
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analysis  problem  in  Section  3 have  two  basic  properties  which  are  extremely 
important  and  useful  for  the  application  of  these  algorithms  to  problems 


^ with  specially  structured  covariance  matrices.  These  properties  are;  I 

!, 

I (1)  the  algorithms  can  take  full  advantage  of  the  special  structures;  and  'i 

I (2)  they  preserve  the  special  structures  throughout  the  solution  process.  |i 

Of  course,  without  either  one  of  the  properties,  the  other  would  not  bo  so  ■ 

useful.  In  order  to  provide  an  example  of  a specially  structured  covariance  * 

matrix,  consider  the  single  index  model  (and  in  general,  a multiple  index 
model  also).  VVe  know  from  (1.5)  that  j 

V = I + Var(Ij^)ej^  j 

which  is  a diagonal  matrix  plus  a rank  one  modification.  This  is  special!  j 

Matrices  of  this  class  (or  even  its  generalization  (1.5))  have  certain  nice  | 

/ i 

properties  which  can  be  used  in  the  new  algorithms  to  facilitate  their  iro-  ^ 

i 

plemontntion.  In  Section  4,  wo  shall  specialize  the  new  algorithms  to  the 
general  index  models  where  we  drop  assumption  (1.4a).  Our  purpose  there  is  ^ 

to  exploit  the  spocia)  .structures  inherent  in  these  index  models  and  to  ^ 

show  how  they  can  bo  used  in  the  algorithms  most  effectively. 


2.  PARAMETRIC  PRINCIPAL  PIVOTING 


2.1.  Background . This  subsection  is  concerned  with  a brief  review  with 
the  background  materials  required  for  the  entire  paper.  For  a given  n- 
vcctor  q and  n X n matrix  M,  the  linear  complementarity  problem,  de- 
noted by  (q/M)  is  that  of  finding  an  n-vector  x satisfying 

(2.1)  q + Mx>^0,  x^O  and  x'^Cq  + Mx)  = 0 . 

Recently,  there  has  been, an  extensive  amount  of  oapers  uublisiied  on  the 
theory  and  applications  of  the  linear  complementarity  problem.  Wo  mention 
three  most  recent  survey  papers  on  this  subject,  Cottle  [3],  Kancko  [14] 
and  Lemko  [18]. 

It  has  been  proved  (see  Samelson  et  nl.  [25])  that  if  M is  a P- 
matrix  i.e.  a real  square  matrix  wit)i  positive  principal  miners,  then  t)ie 
linear  complomontrrity  yiroblem  (q,M)  lias  a unique  solution  for  every  vector 

q • 

The  iiarametrir  linear  comp] ementar ity  proljlem,  denoted  )jy  { (q  + Xp,M)  : 
X c R]  consists  of  the  family  of  linear  complementarity  problems 
(q  + Xp,H)  where  X is  a scalar,  to  he  varied  over  the  whole  real  line, 
and  p is  a given  n-vector.  For  some  introductory  exposition  to  the  pttr- 
ametric  linear  complementarity  problem,  we  refer  to  Cottle  [4]  and  Kaneko 
[15]  . 

Throughout  this  entire  section,  we  assume  that  M is  ci  P-matrix\ 

2 

which  .IS  not  necessarily  symmetric,  and  that  p is  positive  . It  then 
follows  that  the  linear  complementarity  problem  (q  + Xp,M)  has  a unique 
solution  x(X)  for  every  X . In  fact  x(X)  is  a continuous  )3iccewiso  linear 
function  of  X with  finitely  many  breakpoints  (sec  (4]  c.g.).  If  we  let 

(2.2)  T = max  (-q./p.)  , 

l<i<n  ^ 1 

1.  See  also  the  discussion  in  the  last  paragraph  of  this  subsection. 

2.  Of  course,  wo  can  then  assume  for  simplicity  and  without  loss  of  gen- 
erality, that  p is  thc'  vector  of  I's  . However,  we  do  not  asssumo 
this  for  reasons  which  will  become  clear  in  Section  4. 
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then  it  is  clear  that  q + Ap^O  for  all  A^  A.  Thus  x{A)  = 0 for  all 
A A.  To  determine  x(A)  for  those  A less  than  A,  we  resort  to  tlie 
parametric  principal  pivoting  algorithm  which  is  a parametarized  version  of 
Graves'  lexicographic  principal  pivoting  algorithm  [10].  The  parametric 
principal  pivoting  algoritlim  is  valid  if  M is  either  a P-matrix  or  a 
positive  semi-definite  matrix.  It  is  described  in  Cottle  [4] . We  reform- 
ulate the  algorithm  for  the  case  of  a P-matrix  as  follows.  For  those 
readers  who  are  unfamiliar  with  the  theory  of  pivotal  algebra,  we  advise 
them  to  consult  on  the  references  [2],  [5],  [10]. 

Algorithm  1:  Original  version  of  the  parametric  principal  pivoting  al- 

gorithm for  Me  P 

.Step  0.  Determine  A by  (2.2).  Let  the  initial  canonical  tableau  be 
given  by 

u = q + Ap  H-  Mx 


where  q + Ap  > 0 for  all  A>A.  Lot  (u(A),  x(A)) 


(q  -t  Ap,0)  for 


A > A 


Put 


'old 


= X and  go  to  Step  1. 


Step  1 . If  p £ 0 (this  cannot  occur  initially),  let  (u(A),  x(A))  = 
(q  + Ap,0)  on  (-™,  A^^^)  and  terminate.  Otherwise  determine  tlie  new 
critical  value 


now 


ma.x  - : p^  > 0 ) 


and  let  Uie  new  cr it i cal  inde^i  k be  a maximizing  index.  Put  (u(A),x(X)) 

Set  A 


= (q  + Ap,  0)  on  [ A , , , A ] 
^ V/  / Q^d ' new 


old 


A and  to  to  Steo  2 . 

new 


Step  2.  Pivot  on  mj^j,  and  let  u,x,q,p  and  M correspond  to  the  rc'sult- 
ing  canonical  form.  Go  to  Step  1. 

Algorithm  1 may  be  explained  briefly  as  follows.  Initially,  tlio  de- 
termination of  T such  that  q + Ap  > 0 for  all  A^A  is  possible  because  p 
is  )-)Ositive.  This  then  produces  the  solution  x(A)  for  all  A >1  A immedi- 
ately. Next,  one  wants  to  decrease  A and  to  determine  the  solution 
1.  VVe  sometimes  use  tlio  letter  P to  denote  ttie  class  of  P-matrices. 
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x(A)  for  those  A less  than  A . This  is  done  in  precisely 
the  same  way  as  in  parametric  linear  programming,  namely,  by  a ratio  tost. 
In  this  instance,  it  is  a maximum  ratio  test  because  one  is  decreasing  A . 
Having  obtained  the  maximum  ratio  (which  gives  the  largest  extent  for  the 
decrease  of  A , keeping  the  current  basic  solution  satisfy  the  nonnegativ- 
ity conditions)  and  the  maximizing  index,  one  performs  a principal  oivot 
to  replace  the  current  "blocking"  basic  variable  by  its  nonbasic  pair.  One 
then  repeats  the  ratio  tost  to  decrease  A further.  The  algorithm  termi- 
nates when  the  current  p-vector  becomes  non-positive.  In  this  case,  a 
further  decrease  of  X can  result  only  in  an  increase  of  the  values  of  the 
current  basic  varialsles.  Tliorcfore  the  algorithm  terminates. 

Several  remarks  are  important  to  Algorillim  1 and  subsequent  devt lop- 
ment.  Firstly,  the  pivots  in  Step  2 are  alv.’ays  possible  because  of  the 
fact  that  file  matrix  M after  each  such  pivot  remains  a P-matrix  (see  [2] 
c.g.)  thus  has  positive  diagonal  entries.  Secondly,  by  incorporating  a 
lexicographic  rule  to  choose  the  maximizing  index  in  Step  1,  like  the  one 
used  in  Graves'  algorith.m  [10],  the  above  algorithm  will  compute  the  solu- 
tion x(A)  for  all  values  of  A , in  a finite;  number  of  step^.  The  lex- 
icography is  used  to  avoid  degeneracy.  Lastly,  each  pair  of  old  and  nev\' 
critical  values  of  A , nam/'ly,  ^old  ^new'  '’^tually  corresponds  to 

a pair  of  consecutive  breakpoints  of  the  solut  ion  curve  x(A)  . 

As  mentioned  earlier,  the  parametric  principal  pivoting  algorithm 
was  valid  for  a parametric  linear  coTni>lementnrity  problem  with  either  a 
P-matrix  or  a positive  semi -definite  matrix.  In  this  paper,  wc  restrict 
our  consideration  to  the  P-matrix  case  and  leave  the  otlier  case  as  the  r.ub- 
ject  for  anotlier  report. 

2,2.  .Modi  f ica  t ion  and.  extoiis  ion  . Before  proceeding  further,  vjc  explain 
some  notations  to  be  used  later.  Let  M be  an  nx  n matrix  and  let 
I , J c { 1 , . . . n ) . Wo  de  f i no 
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— i 


m . . 


where  I = } and  J = { j , . . . , j ) with  1 < i,  < . . . < i n and 

i-  S It  t s 

1 £ j < . . . < £ n . In  particular,  is  a principal  submatrix  of  M . 

Similarly,  if  q is  an  n-vector,  we  define  q^  = (q.  ) . If 

I = {l,...,n),  we  write  to  mean  Similarly,  if  J = {l,...,n), 

we  write  to  mean  M^j  . Any  vector  indexed  by  the  empty  set  is 

meaiit  to  bo  the  zero  vector.  The  letter  e always  denotes  the  v'ector  of 


In  the  sequel,  we  first  modify  Algorithm  1 in  such  a way  that  it  can 
be  cfficiontiy  adopted  to  large-scale  applications  and  tlien  extend  it  to 
allow  for  upper  bounds  on  tlic  variables.  In  order  to  refrain  us  from  re- 
peating, wo  hereby  declare  tliat  all  the  algorithms  described  in  this  paper 
are  finite  provided  that  degeneracy  or  cycling  is  treated  appropriately  by 
some  lexicographic  r\ilo. 

We  now  e.vplain  how  we  modify  Algorithm  1.  Following  tlie  usual  lan- 
guage in  pivot.al  algebra,  we  call  the  vectors  u and  x in  the  system 


of  linear  equations 


u = q + Xp  + Mx 


Ixi.sjc  and  noii-baslc  vectors  rcnspectivol  y and  their  components  basic  and 
non-bas ic  var iabl es  respectively.  Kow,  suppose  Lliat  .Step  2 of  Algorithm  1 
has  been  performed  several  times  and  that  the  algorithm  is  not  yet  termi- 
nated. Tlie  comjmnents  of  the  current  basic  vector  will  then  consist  of 
some  components  of  tlie  original  x-vector  togetlicr  with  tlie  rest  of  the 
components  of  the  original  u-vcctor.  More  precisely,  tlicrc  are  two  dis- 
joint index  set.s  I and  .T  with  I <■<  j = { l,...,n  ) such  that  I and  J 
consist  of,  resi'cctivoly,  those  indices  of  tlie  components  of  the  original 
u-  and  x-vcctors  that  are  currently  basic.  (Initially,  I = {l,...,n  ) 


I 


and  J - <p  .)  Furtlicrworo,  tho  current  canonical  tableau  can  bo  obtained 
by  performing  a block  pivot  in  tho  initial  tableau.  In  other  words,  writing 
the  initial  tableau  in  partitioned  form  as 


Uj  = qj 

+ Apj 

^II^I 

+ Xpj 

Xj  in 

terms 

of  Xj 

and  u_, 

J 

we  obtain  the 

current  canonical  tableau: 

(2.3a)  Uj  = (q^  - f’j)  + ^ I ~ j”j j‘'’j I ^ 


(2.3b) 

..-1 

- ' 

Letting 

and  hj  be 

t)iG  (unique)  solutions  to  tho 

of  linear  eejuations  respectively, 

(2./!) 

and  defining 

(2.5) 

= - Mjjgj  and  dj  = Pi-Mjjh 

we  may 

write 

(2.3)  as 

(2.6a) 

^I  = 

c + A d J.  + 

(2.6b) 

-9j  - - 

Noting  that  the  ratio  lost  in  Step  1 of  Algoritlim  1 depends  solely  on  thi 
vectors  c, , d,,  g^  and  and  that  performing  the  pivot  in  Stop  2 is 

J.  J.  J J 

actually  equivalent  to  adjusting  tho  sets  I and  J suitably  and  then  to 
obtaining  tho  corresponding  canonical  tableau  (2.6),  we  can  readily  state 
tho  modified  version  of  Algorithm  1 as  follows. 


Algorithm  2:  Modified  version  of  the  parametric  principal  pivoting  for 

M c P . 
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step  0.  Let  X 


Let  1 = ,n}  and  J = C 


Step  1.  Solve  the  systems  of  equations  in  (2.4)  for  and  hj  and 

tlien  compute  and  d^  by  (2.5).  Go  to  Step  2. 

Step  2.  If  and  h >0,  set  x.(X)  = 0 and  x,(X)  = -g_  - Xh. 

X U J.  c/  c 

for  all  X £ ^o2d  terminate.  Otherwise,  determine  the  new  critical 


value 


S)  ■ 

^new  ^ " d~'  i f I n>ax{-  : h . < 0,  j r,  j) } 

i j ^ 

and  let  the  new  critical  index  k be  a maximizing  index.  Put  (X)=  d 

and  x-(X)  = -g  - Xh  for  Xe  IX  ,X  , ,) . Set  X .,  = X and  go  to  Step  3, 
J J new  old  old  new 

Stop  3.  If  k cl,  then  replace  I and  J by  I\(k}  and  ju(k)  re- 

spectively. Go  to  Step  1. 

If  k cJ,  then  replace  I and  J by  I u {k}  and  J \ (k)  I'cspec- 
tively.  Go  to  .Stop  1. 


Rema r k . The  set  J will  always  contain  at  least  one  element  except  initial- 
ly. This  is  clear  from  the  way  the  critical  index  is  chosen  in  Sten  2 and 
the  fact  that  h^  ■ 0 if  J is  a singlt'ton. 

Algoritlim  2 is  basically  the  same  as  Algoritlim  1 except  that  the 
pivot  step  in  the  latter  algorithm  has  been  replaced  by  the  updating  of  the 
two  index  sots  I and  J and  by  solving  the  two  systems  of  linear  equa-,,. 
tions  in  (2.4)  which  always  liave  unique  solutions.  The  ratio  test  is  per- 
formed in  Stop  2 without  any  change.  The  updating  of  the  sots  I and  0 
is  clear  and  requires  no  further  explanation. 

The  idea  of  replacing  the  pivot  steps  by  solving  systems  of  linear 
equation  in  the  above  modification  seems  as  thougli  fairly  simple,  neverthe- 
less it  suggests  a very  efficient  way  to  actually  implement  the  original 
algorithm,  and  is  especially  useful  in  large-scale  applications.  As  a 
matter  of  fact,  the  same  idea  lias  boon  applied  to  Lemke's  algorithm  117) 
for  certain  class  of  linear  complementarity  problems.  The  application  was 
described  in  [26].  It  should  lie  pointed  out  that  the  modification  described 
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in  Alqoritlim  2 is  not  containod  in  [26]  a ] is  obtained  independently  heie. 
Indeed,  tlie  author  is  indebted  to  Profes  Stephen  Robinson  who  has 
called  his  attention  to  the  reference  [2uJ  after  learning  t lie  results  ob- 
tained by  the-  author. 

A large  part  of  the  computational  effort  required  in  Algorithm  2 
consists  of  solving  the  systems  of  linear  equations  in  (2.4)  for  various 
sets  J . Of  course  tht'se  systems  can  be  solved  independently  of  each 
other  by,  for  instance,  factori  :^at  ion  methods.  But  solv'ing  the  systems 
in  tliis  way  is  definitely  not  most  efficient  because-  one  is  not  malting  u.sc 
of  any  old  infoimation  to  generate  new  data.  Here  the  updating  approach 
desciibc'd  in  |26]  can  be-  applied.  The  idea  is  to  avoid  solving  the  whole 
current  systems  of  linear  equations  in  (2.4)  from  scratch,  but  instead, 
by  taking  advantage  of  the  infoimation  one  heir,  about  the?  previous  systems. 
First  note  that  each  time  (2.4)  is  invoked,  the  set  d changes  by  one 
clemonf , or  more  precisely,  increases  or  decroasos  by  the  new  critical  in- 
dex. Tliis  moans  that  the  submatrix  either  increases  or  decreases 

by  one  rov:  and  the  corresponding  column.  If  one  uses  ortlicgonal  (QR) 
factorization  to  solve  tlio  systems  of  equations  in  (2.4)  (see  [11],  e.a.), 

then  liased  on  the  fact  about  the  change  of  the  submatrix  M,,,  it  is  pos- 

u J 

siblo  to  devolojj  fast  numerical  procedures  to  update  the  OR  factors  of  the 

now  M . Tliose  updating  procedures  can  then  be  used  to  facilitate  the 
U J 

cominitation  of  the  solution  vectors  and  li.^  in  (2.4).  The  amount  of 

operations  required  to  compute  t)ie  vectors  g^  and  hj  by  this  approach 

is  of  order  C where  f is  the  size  of  tlic  current  M . Tliis  is  bette 

00 

than  solving  the  systems  of  equations  in  (2.4)  directly  by  an  order  of  mag 
nitude.  The  updating  ))rocodurcs  arc  described  in  more  detail  in  [26].  We 
would  also  like  to  mention  that  updating  formulas  for  the  i^twerses  of 
arc  given  in  [27],  [11].  The  approach  of  using  updated  inverses  to  solve 
the  systems  of  equations  in  (2.4)  is  slightly  inferior  (althougli  of  the 
same  order  of  computational  effort)  tlian  that  of  using  updated  OR  factor 


I 1 


(see  [26]).  Finally,  we  point  out  that  the  total  amount  of  storage  re- 
quired by  using  either  one  of  these  two  approaches  is  of  order  n^  , 

It  is  clear  that  Algorithm  2 has  the  two  properties  mentioned  in  the 
introduction  when  it  is  applied  to  parametric  linear  complementarity  pro- 
blems { (q  + Ap,M) : A e R}  with  specially  structured  matrices  M,  namely, 
it  can  take  advantage  of  and  preserve  the  special  structures  througliout  the 
solution  process.  Sometimes,  certain  special  structure  allows  the  solu- 
tions gj  and  hj  in  (2.4)  to  be  computed  readily  (or  even  explicitly) 
without  the  need  of  going  through  the  updating  procedures  mentioned  in  the 
last  paragraph.  This  is  true  for  example,  in  the  case  of  the  index  models 
as  we  shall  see  in  .Section  4. 

In  the  rest  of  this  subsection,  wo  extend  the  above  analysis  and 
Algorithm  2 to  deal  with  upper  bounds.  A similar  extension  of  Lemko's  al- 
gorithm is  described  in  1261. 

Specifically,  wo  consider  the  following  parametric  linear  comple- 
mentarity problem: 


(2.7) 


u = q + Ap  + Mx  + y ^ 0 


V = a 


> 0 


X 0 
y > 0 


T T 

u X = V y = 0 


where  A,q,p  and  M are  as  above  and  a is  a given  positive  vector  of  upper 

bounds  on  the  varia)3les.  Problems  of  the  tyoc  (2.7)  arise,  for  example, 

from  quadratic  minimization  oroblems  with  unper  bounds.  For  fixed  value 

of  A,  problem  (2.7)  is  a linear  complementarity  problem  v;ith  vector 

^q  + Apj  matrix  (_^  q)  . Tliercfore,  by  varying  A over  the  whole 

a 

real  line,  problem  (2.7)  becomes  the  parametric  linear  complementarity 


problem  { C^)  + q' ^ hinear  complementarity  problems  of 


the  tvpe  (2.7)  with  fixed  A have  been  studied  intensively  in  two  recent 
papers  16],  I2l] . See  also  (16).  In  particular,  it  was  proved  in  (6]  that 
if  M is  a P-matrix,  then  for  every  fixed  A,  problem  (2.7)  has  a solution 
x^(\)  (from  which  the  corresponding  y'NX)  can  bo  generated  easily).  The 


/ I 
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uniqueness  of  x®(,\)  was  oroved  in  121].  Our  obiective  here  is  to  extend 
Algorithm  2 to  determine  x ®(X)  as  a function  of  X parametrically.  For  a 
unified  study  on  problems  (2.7)  with  X = 0 and  (2.1),  we  refer  to  122] 
where  the  author  treats  the  vector  a of  upper  bounds  as  having  some  com- 
ponents equal  to  infinity. 

The  reason  that  Algorithm  2 is  not  directly  applicable  to  problem  (2.7) 
is  o)3vious:  (^^  J)  is  not  a P-matrix  even  if  M is.  Suppose  that  we 

have  obtained  some  value  X of  X for  which  we  know  the  solutions  x^(X)  for 
all  X>^  X and  now  we  want  to  decrease  X to  complete  the  solution  process. 

In  the  solution  x^(X),  some  components  will  be  basic  and  others  nonbasic. 
Among  tlie  basic  components,  some  will  at  their  upper  bounds.  Let  I and 
J denote  respectively,  the  sets  of  indices  of  the  now  basic  u-  and  x- 
variables,  with  J consisting  of  and  where  corresponds  to 

tlie  subset  of  J such  tliat  the  basic  x-variables  with  indices  in  J2  are 
at  tliei.r  upper  bounds  and  is  tlic  complement  of  in  J.  (Initially, 

X is  given  liy  (2.2),  I = (I,...,  11}  , and  = J2  = <f' . ) With  the  sets 

I,  and  J2,  we  caj\  v.’rite  tlie  initial  canonical  tableau  in  partitioned 
form  as: 


“l 

= ‘Ij  + XPj  H 

"ll^l  + ^ ^ 

= q-T  + xp. 

V, 

= a,. 

” >^T 

I 

3 

I 

V 

= ci 

“X  , 

'^1 

■^1 

+ y- 


and  the  intermediate  canonical  tableau  before  tlic  furtlier  decrease  of  X as: 


-15- 


1 


"M  ^ V 

jij/ai 

Xy  = a -V, 

*^2  *^2  ■^2 

Vj  = 3y  -Xj 

*~M^  M V +M^v 

JiJi 

+ M m”^  V 

J2J, 

Letting  g and  h.,  denote  the  (unique)  solutions  to  the  following 

Jl  Jj 

systems  of  lineai-  equations  respectively, 

(2.8)  M g =q.,  + M a and  M h -D, 

Jl^l  JjJ2  J2 

and  defining 

~ ■*■  ’'’lJ2^J2'  ^1  '’l  " 


(2.9b) 


c,  = q, 

•^2  *^2 


”■^2^1^'^!  ^ ”‘'2''2^-^2'  ‘^‘'2  “ ^*^2  ■ "'^2'^1^'^1 


we  may  rewrite  the  tableau  corresponding  to  X - T as 
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(2.10) 


Uj  = Cj  + Adj 


" "9j  ■ 

Jl 

X = a 
^2  ^2 


Vj  = aj 


V = (a,  + g_  ) + Xh^ 

Jl  Jj 

y,  = -c  - Ad 

■^2  ^2  ^2 


+ nonbasic  portion. 


It  is  now  clear  from  (2.10)  how  to  proceed  to  determine  the  decrease 
of  A . Based  on  (2.30)  and  the  same  idea  as  before,  wo  deduce  the  next 
algorithm  which  is  an  extension  of  Algorithm  2 to  solve  problem  (2.7). 

Algorithm  3;  Extension  of  the  parametric  pi'incipal  pivoting  algorithm  to 
deal  with  upper  bounds  for  M e h . 


Step  0 . I.et  A 


Let  I = {l,...,n)  and  = J2  = <^  . 


Step  1.  Solve  the  systems  of  equations  in  (2.8)  for  q,  and  h,  and 

- Jl  Jj 

then  compute  c_,  d , c and  d by  (2.9).  Go  to  Steo  2. 

1 1 J 2 *^2 

step  2.  If  dj£0,  hj  =0  and  dj  > 0,  set  x®(A)  = a for  all  A £ ^old 
and  terminate  (in  fact,  I and  will  then  bo  empty).  Otherwise  determine 

the  new  critical  value 

c . 

A„„  , = raax{  max  {-  d . > 0,  i c J } , max  (-  c : h.  <0,  j,EJ,)  , 

new  d.i  h.  i,  11 

1 3^  -'1 


max  {“ 


a . +q . 

3l  3l 


: h.  max{-  3 : d.  <0,  j2t0r^)l 

3i  h \ 


and  let  the  new  critical  index  )t  be  a maximizing  index  . Put  x“j(A)  = 0 , 

(A)  - -g.,  - Ah  and  x®  (A)  = a.,  for  Ac  IA^,.,A  1 . Set 

»J  ^ J J 0 2^  oXd  now 

A , , = A and  go  to  Step  3. 

old  now 

Step  3.  If  k e I,  then  replace  I and  J^  by  I X (k)  and  Jj  t>  {);) 
respectively.  Go  to  Step  1. 

If  kf  J2,  then  replace  and  J2  by  Jj^o{k)  and  J2X  (k) 

re.speet ively . Go  to  step  1.  -17- 


If 


by 


I u {k}  and 


2 by  \ {k ) and 


k c and  hj^  < 0,  then  replace  I and 
V {k}  respectively.  Go  to  Step  1. 

If  k e and  hj^  > 0,  then  replace  and  J 

J2'-’(k)  respectively.  Go  to  Step  1. 

The  termination  of  Algorithm  3 deserves  some  more  discussion.  Ac- 
cording to  Step  2 of  the  algorithm,  wo  have  at  termination  d < 0,  h =0 

1 - 


and 


>^2- 


> 0,  and  for  all  X < X 


- "old' 


all  the  components  of  the  solution  x (X) 


are  at  their  upper  bounds.  In  fact,  h. 

Therefore  since  p is  positive,  we  must  have  I = . On  the  otlier  hand,  the 


...  =0  if  and  only  if  J.,  = 

1 ^ 


fact  that  the  solvition  x®(\)  = a for  all  X ^old  expected  and  can 

be  justified  directly  from  the  original  formulation  (2.7)  of  the  problem. 

As  we  know,  the  solution  x(X)  is  a piecewise  linear  function  of  X 
Therefore,  there  exists  a value  X (which  is  in  fact,  X , , at  termination 
of  Algorithm  such  that  x^(X)  is  linear  for  all  X < X.  Now  since  x^(X) 
is  bounded  above  and  below,  it  must  bo  independent  of  X for  X < X,»  i-e-  x^(X) 
is  equal  to  some  constant  vector  for  X ^ X^.  (This  is  not  true  in  the  case 
witliout  tlic  upper  bounds.)  Suppose  the  i-th  component  of  this  vector  lies 
strictly  below  its  upper  liound . Then  we  must  have,  by  complementarity 


0 £ = qp  + Xp^  + (Mx  (X))^ 


which  certainly  can  not  bo  true  for  all  negative  X . Therefore  it  follows 
that  for  all  X _<  X , x^(X)  is  equal  to  the  vector  a of  upper  bounds. 

It  should  be  pointed  out  that  every  time  Step  3 is  passed,  tlie  set 
changes  (by  one  clement)  . It  may  sometimes  be  empty  (as  it  must  bo  at 
termination).  Kurt  liermoro , the  updating  procedures  described  earlier  can  be 
used  to  solve  the  systems  of  equations  in  (2.8)  more  efficiently. 

We  conclude  this  subsection  and  also  this  section  by  remarking  that 
under  suitable  modification.  Algorithm  3 can  bo  applied  to  the  situation 
where  some  of  the  upper  bounds  arc  infinity.  Roughly  speaking,  this  is  done 
by  restricting  the  set  J2  to  include  only  those  indices  for  which  the 
upper  bounds  are  less  than  infinity.  See  the  discussion  in  122]  for  the 


3 . APPLICATIONS  TO  PORTFOLIO  ANALYSIS 

3.1.  The  case  without  explicit  upper  bounds.  In  this  subsection,  we 
apply  Algorithm  2 to  develop  a new  and  efficient  algorithm  (Algorithm  5)  for 
the  following  parametric  quadratic  program. 


(3.1) 


. It  t 

minimize  ^ x Mx  + (q  + Or)  x 


subject  to 
and 


X ^ 0 

T 1 
p X = 1 


where  M is  a given  symmetric  positive  definite  nXn  matrix;  p,q  and  r 
are  given  n-vectors  with  p positive;  and  0 is  a nonnegative  parameter.  We 
recall  that  the  general  portfolio  analysis  problem  with  a nonsingular  co- 
variance  matrix  is  a problem  of  this  type.  Furthermore  v.'e  have  q = 0 and 
p = e (see  (1.2)).  The  first  identification  (q  = 0)  has  certain  ad- 
vantage in  the  implementation  of  the  first  algorithm  proposed  below. 


Since  p is  positive,  the  variables  in  problem  (3.1)  are  always  im- 
plicitly bounded  above.  Therefore  since  M is  symmetric  and  positive  def- 
inite, problem  (3.1)  has  a unique  solution  x(P)  for  every  0 . Further- 
more, x(6)  is  also  the  unique  vector  satisfying  the  following  Kuhn-Tuc)'.er 

conditions 

(3.2a)  u = q + Ap  + Or  + Mx  _>  0 , x^O 

(3.2b)  u"^x  = 0 

(3 . 2c)  p"^x  = 1 

where  x is  the  Lagrange  multiplier  associated  with  the  equality  constraint 
(3.2c)  and  is  dependent  on  6 . Also  X is  not  restricted  in  sign.  Consider 
0=0  first.  The  resulting  problem  is 


(3.3a) 

U = q + Xp  + Mx>0, 

X > 0 

(3.3b) 

T 

U X = 0 

(3.3c) 

T , 

p X = 1 . 

Now  problem  (3.3)  can  be  regarded  as  the  parametric  linear  comple 
mentarity  problem  { (q  + Xp,M):  X r R)  coupled  with  the  single  equality 
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constraint  (3.3c).  In  other  words,  to  obtain  x(0),  we  can  let  X be 


a 


I 

t 

1 


I 

parameter  varying  over  the  real  line  and  aim  at  deteriniiiincj  a suitable 

* * 

A such  that  the  solution  x(A  ) to  the  linear  complementarity  problem 
(q  + X p,M)  satisfies  condition  (3.3c)  also.  It  is  then  clear  that  x(0) 

IS  given  precisely  by  x(X  ) . One  way  to  find  such  a X and  the  solution 

* 

x(X  ) is  to  solve  the  parametric  linear  complementarity  problem 
{ (q  + Xp,M):X  c r)  and  to  develop  a systematic  way  to  check  if  the  solution 
x(X)  satisfies  condition  (3.3c).  Algorithm  2 can  be  applied  here  to  obtain 
x(X)  for  every  X and  in  fact,  it  also  provides  an  easy  way  to  check  (3.3c). 

In  what  follows,  we  explain  how  this  latter  part  can  be  done.  Roughly  speak- 
ing, it  is  achieved  by  linear  interpolation. 

Recall  that  x(X)  is  a continuous  piecewise  linear  function  of  X 
Indeed,  according  to  Step  2 in  Algorithm  2,  we  have 

X (X)  = 0 and  x,(X)  = -g  - Xh  for  Xr  [X  ,X  ] 

I J ■’J  J new'  old 

where  and  X^^^^  correspond  to  a pair  of  consecutive  breakpoints  of 

the  solution  curve  x(X);  and  h.,  arc  solutions  to  the  systems  of  linear 

U J 

equations  in  (2.4).  Thus, 

p’’x'(X)  -pj  g^  - XpJ  hj  , 

T— 

for  X f |X  ,X  ,,).  Setting  p x(X)  = 1 gives 
now  o J.  c. 

(3.4)  1 = -(1  + pj  gj)/pJhj  , 

j'rovided  that  .T  is  nonempty.  Geometrically,  the  quantity  "X  given  by 

(3.4)  is  precisely  tliat  value  of  X wliere  the  (riossibly  extended)  line  seq- 

ment  of  t Ik  solution  curve  x(X)  in  the  interval  IX  ,X  , ,)  meets  the  line 

new  old 

p x(*)  - 1 . So  i f this  quantity  "X  lies  within  the  interval  (i.e.  if  the 

* 

intcisection  occvirs  within  the  interval),  we  have  found  X = X . Other- 
wise we  decrease  X and  proceed  to  the  next  interval.  Eventually  in  a finite 

* 

numbei  of  steps,  this  process  will  terminate  with  a X and  of  course,  the 
snluticin  x(X  ) (=  x(0))  . We  summarize  the  discussion  in  the  next  algorithm 

which  solves  problem  (3.1)  witli  0 = 0. 
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Algorithm  4;  QP  with  0=0  and  no  upper  bounds. 


Step  0. 

Determine 

X - , = max  {-  — } 
l<i<n  Pi 

and  let 

k be  a 

maximizing  index 

Lot  I = 

(1, . . . ,n) 

X (k)  and  J = {k} . 

Go  to  Step  1 . 

Step  1. 

Solve  the 

systems  of  linear  equations 

in  (2.4) 

for  g.,  and  h. 

and  then 

compute  c 

J and  dj  by  (2.5) 

. Go  to 

Step  2. 

Step  2. 

Compute 

X by  (3.4).  If  dj 

< 0 and 

then  set 

(3.5) 

X*  = X,  Xj (0)  = 0 

and  X 

ij 

(0)  = -g^ 

- X*hj 

and  terminate.  Otherwise  determine 

c . g . 

(3.6)  1 = max  {max  {-  : d,  > 0,  i c I),  max{-  ~:h.<0,  jcJ)} 

new  u . 1 n . 1 

1 n 

— * /V 

If  Xc  J ^riew' ^old^  ' define  X and  x(0)  by  (3.5)  and  terminate. 

Otherwisr  let  k be  a maximizing  index  in  (3.6).  Go  to  Step  3. 

Step  3.  If  k r I,  replace  I and  J by  I \ (k)  and  JU  {k}  respective! 
Go  to  Ste]5  1. 

If  k r,  J,  replace!  I and  J by  I U {k}  and  J \ (k)  respectively. 
Go  to  r-tep  1. 

Hemark  1 . As  in  Algorithm  2,  the  set  J always  contains  at  least  one  ele- 
ment, so  that  the  guantily  X in  (3.4)  is  always  well-defined  by  the  posi- 
tive definitene.sG  of  M . 


Remark  2.  If  q happens  to  be  zero,  as  in  t))o  portfolio  analysis  problem 

(1.2),  then  it  is  easy  to  see  that  X , , and  X will  always  be  zero. 

old  new 

Therefore  t)je  ratio  test  (3.6)  and  the  one  in  Stop  0 can  be  eliminated. 
Furtliermore , the  maximi  ziiig  index  k can  be  cho.sen  "arbitrarily"  (subject 
to  non-cycling)  provided  cither  one  of  the  following  two  conditions  is 
satisfied : 


(1) 

(2) 


dj^  > 0,  k c I , 

)ij^  < 0,  k c J 


In  Step  0,  k can  be  chosen  arbitrarily  in  {l,...,n)  . Finally,  if  q = 0, 
Algorithm  4 terminates  if  and  only  if  and  ^>j  some  suitable 


I and  J . 
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Remark  3 . The  assumption  that  M is  positive  definite  is  essential  for 
the  success  of  Algorithm  4 to  produce  the  solution  x(0)  to  problem  (3.1) 
with  0=0.  This  is  because,  on  the  one  hand,  the  assumption  guarantees 
the  existence  and  uniqueness  of  x(0)  which  also  satisfies  conditions 
(3. 2a) - (3. 2c)  for  0=0  and  suitable  1;  on  the  other  hand,  it  guarantees 
that  for  every  X,  the  linear  complementarity  problem  (q  + Xp,M)  has  a 
unique  solution.  Consequently,  by  solving  the  parametric  linear  comple- 
mentarity problem  { (q  + Xp,M) : X e R},  one  is  bound  to  arrive  at  the  X* 
and  the  solution  x(0)  given  in  (3.5). 

Using  Algorithm  4 as  a start,  we  may  proceed  to  develop  a parametric 
algorithm  for  problem  (3.1),  or  equivalently,  problem  (3.2).  Observe  that 
there  are  "two"  parameters  X and  0 in  (3.2),  with  X depending  0 . It  will 
be  shown  in  wliat  follows  that  X can  alv;ays  be  eliminated  by  means  of  the 
equality  constraint  (3.2c)  and  exprof^scd  in  terms  of  the  values  of  the 
current  basic  x-variables  and  the  real  parameter  0 . Hence  we  are  left  with 
one  parameter  and  can  therefore  proceed  as  before. 

Let  I and  J be  the  index  sets  at  termination  of  Algorithm  4. 

Then  we  can  write 


u = q + Xp  + Or  + Mx 

in  canonical  form  v/ith  respect  to  I and  J as: 


(3.7a)  Uj  = Cj  t Xdj  t Obj  t (Mjj  - Mj^M-’m^^)Xj  4 
(3.7b)  Xj  = -g_^  - Xhj  - Ofj  - 4 


where  f , g and  li  are  (unique)  solutions  to  the  systems  of  linear  equa- 

J u u 

tions 


(3.8) 


respectively;  c^  and  d^  are  given  in  (2.5);  and  finally 
Setting  the  nonbasic  variables  x^  and  Uj 


at  zero,  we  liave 


-Pj9j  ' - ®Pjfj 


Setting  p x = 1 gives 


T T 

l+P-,g-r  ptf, 

(3.9)  X = - 0-4^  , 

Pj^j  Pj^j 

provided  that  J is  nonempty.  Substituting  X into  (3.7),  we  obtain 


(3.10a) 


(3.10b) 


T T 

i+p-,g-,  p-jf. 

(Cj  - — + 0(bj  - -^d^) 

Pj'b  Pj’'a 


i+Pib-,  Pjfj 

(yj  - - 4-  hj) 

Pj'\t  Pj^J  . 


+ nonbasic  portion 


Starting  with  0 = 0,  wo  may  use  tableau  (3.10)  to  determine  the  in- 
crease of  0 and  proceed  as  before.  Recalling  the  definition  of  Cj , dj. 
and  bj. , wo  may  write 

Cj  - dj  (q^  - ijr  Pj)  - Mj.j(gj  - hj) 

Pj^^J  Pj”j  Pj'^J 

T,,  T,.  T, 

, Pj^J  . _ , . _ Pj^J  . M /<=  _ V.  ^ 

“ T T,“  Pl^  T. 

Pj^J  Pj'b  Pj^^J 


Now  we  can  formulate  the  parametric  algorithm  for  problem  (3.1)  as 
follows . 

Algorithm  5:  Parametric  QP  with  no  upper  bounds. 

Step  0.  Use  Algorithm  4 to  obtain  x(0)  and  let  I and  J be  the  two 
index  secs  at  termination  of  the  algorithm.  Let  0^]^^  = 0 • 

Step  1 . Solve  the  .systems  of  linear  equations  in  (3.8)  for  f..,  g . and 


hj  . Compute 


t 


(3.11) 


rp  m 


X = - (1  i Pj9j)/Pj\y  = - 


(3.12) 

= 

+ vhj 

(3.13) 

Cj  = (qj  + Xpj)  - 

”lJ®J' 

dj  = (rj 

+ npj)  - Mjjtj  . 

Step  2. 

If 

d^  > 0 and  t < 0 , 

set 

Xj(0)  = 0 

and  x-(0)  -s.,  for 

a J 

— old 

and 

terminate  (in  fact, 

will  then 

be  zero) . Otherwise  de- 

termine 

(3.14) 

0 

new 

c . 

= min{min{“  d, 

d . 1 
1 

< 0,  i 

El),  min { 

s . 

- -1;  t , > 0,  j E J}  ) 

j ^ 

and  let 

k be  a minimizing  index.  Put  Xj.(0)  = 

0 and  X (B)  = -s.,  - Ot, 
u i3  a 

for  0 c 

'«old 

, 0 ] . Sot  0 , , = 

new  old 

0 

now 

and  go  to 

step  3. 

Stop  3. 

If 

k e I , replace  I 

and 

J by  I Wk)  and  Jt’lk)  respectiv 

oly.  Go 

to  step  1. 

If  )c  cj,  replace-  I and  J by  I >■'  (k)  and  .7  \{k)  respectively. 
Go  to  Step  1. 

Komar k 1.  The  se-t  J will  never  be  empty  because  of  the  fact  that  t^  = 0 

for  singleton  J and  also  of  the  way  the  minimizing  index  k is  chosen. 
Therefore,  the  quantities  X and  p in  (3.11)  are  always  well-defined. 

Remark  2 . This  time,  <'i  minimum  ratio  test  is  performed  in  Stop  2 because 
we  are  increasing  0 . 

Remark  3.  The  following  two  identities  are  clear  from  (3.11)  and  (3.12) 
(3.15)  97)d  P^t  T = 0 

J J J 

They  guarantee  tliat  for  every  0,  the  solution  x(0)  defined  in  Stop  2 
T /N 

satisfies  p x(0)  =1  . 

Remark  4.  As  different  from  the  last  algorithm,  the  condition  (q  = 0) 
does  not  induce  too  much  simpl ication  in  Algorithm  5 except  possibly  in  the 
compu  t a t i on  s ( 3 . 1 1 ) - ( 3 . 1 3 ) . 
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It  is  worth  point iiuj  out  two  rather  interesting  consequences  at  the 
termination  of  Algorithm  5.  The  first  consequence  is  that  the  solution 
x(0)  for  all  G greater  than  or  equal  to  the  value  of  0 at  termination  is 
equal  to  some  constant  vector  which  is  independent  of  the  parameter  0 . The 
direct  justification  for  this  is  similar  to  that  for  Algorithm  3 which  we 
have  discussed  earlier.  The  second  consequence  reflects  certain  rolationshi 
betweeji  the  solution  x(G)  at  termination  and  the  ratios 

cording  to  Step  2,  we  have  d > 0 and  t,<  0 at  termination.  Since  p is 
positive,  (3.15)  then  implies  that  indeed  t = 0 . Recalling  the  defini- 
tions  of  and  tj,  we  may  deduce  tluit 

(3.1G)  r.^  t )ip  > 0 and  i-  + pp  = 0 

or  equ.i vnlently , 


X . j-  . 

(3.16) '  min{--i)  for  every  jcj  . 

TJiis  latter  relation  .says  that  if  0 is  large  enougli,  then  component  j of 
the  solution  x(G)  can  not  be  i^ositive  unless  the  corresponding  ratio 
I'j/Pj  is  no  greater  tlian  al  vhe  other  ratios  (r^/p^)  . 

3.2.  The  case  with  explicit  upper  bounds.  In  tliis  subsection,  we  extend 
the  analysis  and  algoritlmis  in  tlic  previous  subsection  to  tlie  problem  with 
explicit  upper  bounds; 

IT  T 

(3.17)  minimize  ^ x Mx  + (q  1 Or)  x 

.subject  to  a > X > 0 

T 

and  p X = 1 


where  p,q,r  and  M are  ns  above  and  a is  a positive  vector 
a^  < 1/p^  for  every  i . This  latter  condition  is  imposed  on 

so  that  the  upper  bounds  will  not  be  redundant.  Furthermore, 

T 

the  problem  (3.17)  to  be  feasible,  we  must  have  P a^l  . If 
then  every  feasible  vector  x must  satisfy  x - a because 


sucli  that 
the  vector  a 
in  order  for 

T , 

P a = 1, 

p is  positive, 
p"^a  > 1 


The  problem  (3.17)  thus  becomes  Irivial.  Consequently,  we  assume 


M 


throughout  the  discussion  below.  The  portfolio  analysis  problem  with  uoper 
bounds  (1.7)  is  of  the  form  (3.17)  with  q = 0 and  p = e . 

For  every  fixed  0,  let  x^(0)  denote  the  (unique)  solution  to 
problem  (3.17).  Since  M is  symmetric  and  j-iositivo  definite,  problem  (3.17) 
is  actually  equivalent  to  its  Kulin-Tucker  optimality  conditions 


(3.18a) 
(3.18b) 
(3.18c) 
(3. 18d) 


u=q+Xp+Or+Mx+y  ^0, 

V = a - X i 0/ 

0 


X > 0 
y >0 


T 

u X 


T 

V y 


T 1 

p X = 1 


Consider  0 - 0 . Then  condition  (3. 18n) - (3. 18c)  can  be  regarded  as 
a parametric  linear  complementarity  problem  of  the  form  (2.1)  with  >,  as  the 
parameter.  Therefore:,  in  order  to  comi^uto  x^(0)  wliich  must  satisfy  con- 
ditions ( 3 . 1 8:i)  - (3 . 18d)  fcir  some  suitable  X , we  may  use  Algorithm  3 and 
incorporate  a rule  to  checl;  condition  (3.18d).  Tlio  idea  is  exactly  the 
same  for  the  case  with  no  explicit  upper  bound,-.;  therefore  we  omit  the  de- 
tail. Tlic  next  algoriihiii  extends  Algoritiim  4 to  solve  problem  (3.17)  with 
0 =0  . 


Algoritiim  6:  QP  wit!  0 


and  exjilieit  upper  1,'ounds. 


Detc'rminc  be  a maximizing  index. 

1 ^ i < n ^ 


Step  0 . 

Let  1 = {],...  ,n)  \ {k  },  = ik  ) and  ~ ‘I  ■ Go  to  .8tc)-''  1 . 


Step  1 


Solve  the  systems  of  linear  equations  in  (2.8)  for 


and 


and  then  compute  c, , d , c and  d,,  by  (2.9).  Go  to  Step  2. 
J-  J-  ^ ^ ? 


Step  2, 
(3.19) 


Determine 


c . 


'’ll 


= max  (max  {-  : d.  > 0,  i ( 1 ),  max  {-  v- — - : h.  <0,  j,  ^ ^ t 

d.  j ' h.  j,  -'1  1 

1 -’1 


a . +g  . 
3l 

max  {-  — -7- 


c . 
32 


3l  , r 

--j- : h.  > ),  maxi-  ^ ^ ^ 

3l  32  2 


d 


0,  j c j )} 


Case  (i)  *3^^  t (Ji . Put 


I 


X = -(1  + Pj  q - Pj  a )/p^  h . 

Jl  J2  J2 

If  X c [X  fX  , , set  X = X , 

new'  old 

x®(0)  = 0,  (0)  = -g^  - X h,  and  x^  (0)  = 

I J2  J2 

and  terminate.  Otherwise  let  k be  a maximizing  index  in  (3.19).  Set 

X_.  , = X and  go  to  Stop  3. 

old  now 

T 

Case  (11)  J-  = (fi  (or  equivalently,  h_  =0).  If  p_  a.  = 1,  set 

'^l  “^2  '^2 

X = X 

now 

x^(0)  = 0 and  x^  (0)  = a 
^ ^2  ^2 

and  terminate.  Otlierwise,  let  k be  a maximizing  index  in  (3.19).  Set 

X , , = X and  go  to  Step  3. 

old  new  ^ ^ 

Step  3.  If  k r,  I,  replace  1 and  by  I ' {k)  and  *-)  {)■.}  respec- 

tively. Go  to  Step  1. 

If  k c J2,  replace  and  d2  '^2  ^ 

spcctivoly.  Co  to  Step  1. 

If  k t dj  and  hj^  < 0,  replace  I and  by  I 'J(k)  and 

\ {k}  respectively.  Go  to  Step  1. 

If  k c and  )ij^  > 0,  replace  Jj  and  J2  by  {);)  and 

J2tt{k)  respectively.  Go  to  Step  1. 

The  oxi.stcncc  of  a maximizing  index  k in  (3.19)  still  requires 

justification.  However,  ttiis  is  clear  because  if  d,.  < 0,  h-  =0  and 

I - 

d >0,  then  according  to  wliat  we  liave  remarked  earlier,  those  conditions 
indicate  that  in  fact,  1 - = <(i  and  that  we  have  arrived  at  tlie  last 

breakpoint  of  the  solution  curve  x®(X)  . These  latter  two  facts  together 

would  imply  x^(0)  = a,  contradicting  p'^'a  1 1 . Therefore  a maximizing 
index  k in  (3.19)  always  exists. 

Finally,  we  present  tlie  parametric  algorithm  to  solve  pinblo.m  (3.17) 
It  is  an  extension  of  Algorithm  5 to  deal  with  explicit  upper  bounds. 
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Algorithm  7:  Parametric  QP  with  explicit  upper  bounds. 

Step  0.  Use  Algorithm  6 to  obtain  x®(0)  and  let  I,  and  J2  bo  the 

three  index  sets  at  termination  of  the  algoritlim.  Let  6 = 0 . If 

old 

_ * _ 

= <t'  « set  X = \ , p = 0 and  go  directly  to  (3.23)  . Otherwise  go  to 


Step  1. 

Stop  1. 

Solve  the 

(3.20) 

”j^jj 

for  f , 
•^1 

^7 

“^1 

and  h 

(3.21) 

I 

= - (1 

(3.22) 

^^1 

= 07  + 

(3.23) 

"•l 

= (qj  + 

(3.24) 

•^2 

= (qj 

Step  2. 

If 

d^  > 0, 

l‘"2  ^2 


'1  "1 


2 2 


1 ^1 


’1  '^1  1 1 


Jl' 


2 2 1 1 


S-,  , d^  = (Sj  + ppj) 


IJl  .7^ 


"2"  2 ' 2 2"1  •-1 


J2J1  J,' 


and  X (0)  = a,  for  0 >0  - , and  terminate.  Otherwise  determine 
'J  2 2 ~~  o 1 ci 


(3 


^ 1 

.25)  0^^^  = min  {n.in  {-  d . < 0 , del),  mi)'  {-  ^ — : t.  >0,  < 

i jl  ^3 

a . 7 s . c . 

=>1  3,  , , 3o 

■r t.  <0,  HiC  min{-  : d.  >0,  j,cJ, 

i -*1  ■*  “-1^2  ^ 

3l  A 32 


min{- 


)} 


and 

let 

1:  be  a minimizing  index.  lUit 

Xi( 

and 

^a 

^7 

*^2 

(0)  - for  Ot 

Set 

Step 

_3. 

If  c I»  replace  I and  Jj^ 

by 

ly. 

Go 

to  Step  1 . 

= 0,  (e) 

'^1 

n , j 0 
old  new 


= -s  - at 

Jl 

and  go  to  Step  3. 
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1 


I 


If  k c ^2'  replace  and  J2  f*y  and  J2  \ {k)  respective- 

ly. Go  to  Stop  1. 

If  k c Jj^  and  tj^>0,  replace  1 and  by  iu{k}  and  (k) 

respectively.  Go  to  Step  1. 

If  k t and  < 0 , replace  and  J2  by  \ {k}  and  J2U{k} 

respectively . Go  to  Step  1. 

Remark  1.  The  set  0^^  will  never  be  empty  except  possibly  at  the  termina- 
tion of  Algorithm  C.  Tlie  reason  is  the  same  as  in  the  case  with  no  upper 
bounds.  Kurthermure  (3.15)  holds;  if  J is  replaced  by  . 

Remark  2.  The  following  analog  of  relation  (3.16)  holds  at  termination 


of  Algorithm  7 : 


(3.26) 


or  cquivalc-nt ly , 


r . r . 

- ^1  ^2 

(3.26)  ' min{  — -1  ^ " V‘  ^ f ^ for  every  j.^  cd, 

ir.l  I’i  Pjj  .i2‘^*2  ^’^2 


This  Inttcr  rolat.ion  (3.26)'  indicates  liow  the  so]ulion  (O  ) is  related 

to  1 lie  ratios  (r./f’*)^  for  su f f icic'nt.1  y larqo  0 . 

^ ^ i-l 


4. 


SPECIAI.I7.ATI0NS  TO  INDEX  MODELS 


t 

i 

I 

I 


i 


In  this  section,  we  apply  the  algorithms  developed  in  the  last  sec- 
tion to  some  specific  models  in  portfolio  analysis,  namely,  the  index  models. 

Recall  that  in  an  m-index  model,  the  returns  of  the  securities  are  given  by 
(1.3)  and  there  are  three  assumptions  ( 1 . 4a) - (1 . 4c)  concerning  various  co- 
variances. Since  it  is  rather  unreasonable  in  practice  to  assume  that  no 
two  indices  are  correlated,  wo  will  drop  the  assumption  (1.4a)  throughout 
tlie  following  discussion.  We  denote  by  C the  covariance  matrix  of  the 

indices,  i.e.  C = (c. .)  ^ whore 

1 3 ni  X m 

c^j  = cov(T^,lj)  for  every  i,j  = l,...,m 

Under  tliis  setting,  it  is  easy  to  deduce  that  the  covariance  matrix 
V of  portfolio  retains  is  given  by 

(4.1)  V = >;  + (pj^,  . . . ,Bj^)  C(p^ 

where  K is  the  nonnegative  diagonal  mati'ix  whose  i-th  diagonal  element  is 

equal  to  Var(c^)  and  p.^  = (p^^  j = l,...,m  . If  C is 

diagonal,  then  (4.1)  reduces  to  (1.5)  mentioned  in  the  introduction.  Note 

that  (P,,...,P  ) is  an  nXm  matrix. 

1 m 

In  the  sequel,  we  add  the  following  assumption  to  the  model: 

(A  ) the  diavgonal  entries  of  ''  are  all  positive,  i.c.  Var(i^)  >0  for 

i = 1 , . . . , n . 

Since  the  can  be  viewed  as  the  differences  betv.’ecn  the  actual 

security  returns  and  those  obtained  by  taking  linear  combinations  of  the  in-  1 

dices,  the  assumption  (A)  says  that  tliere  are  always  errors  in  measurinq  the 
actual  security  returns  by  linear  combinations  of  the  indices.  Fortunately, 
this  is  not  a too  unreasonable  assumjition  because  the  coefficients  of  com- 
binations, i.e.  the  are  usually  obtained  by  statistical  estimation 

based  on  past  lii story  in  the  stock  market  , thus  are  themselves  subject  to 
errors.  Presumably,  taking  linear  comljinations  involving  those  erroneous 
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B^j's  would  provide  only  approximate  and  not  exact  values  of  the  actual 
security  returns.  Unfortunately,  assumption  (a)  will  exclude  ri.skless 
securities  in  the  model. 

If  the  index  models  are  of  any  significance  in  reducing  the  compu- 
tational and  informational  complexity  of  the  general  portfolio  analysis 
problem,  it  would  seem  very  reasonable  not  to  include  too  many  indices  in 
the  models;  or  in  other  words,  we  may  consider  m relatively  small  com- 
pared with  n . This  latter  condition  is  particularly  meaningful  when  n 
is  very  large. 

Under  assumption  (a)  , the  matrix  V given  in  (4.1)  is  symmetric 
and  positive  definite.  Thei'efore  the  algorithms  developed  in  Section  3 are 
applicable.  Of  course,  (A)  is  only  a sufficient  condition  for  V to  bo 
positive  definite.  Making  the  change  of  variables, 

i 

2 = f X 

where  J,  is  the  diagonal  matrix  whoso  dingona.1  entries  are  the  sgucire  roots 
of  t)ie  corresponding  diagonaJ  entries  of  V.  , we  In^iy  reformulate  problem 
(1.2)  ((1.7))  as 


(4.2) 

minimi 2c 

It  t 

2’  z Mz  4 Or  * 2 

subject  to 

(a  ' 2l)  z > 0 

and 

p'"*-/.  - 1 

where 

(4.3) 

M = I + (B',..., 

-i 

6,:/ 

(4.4a) 

P!  - Bi 

i 

i = 1 , . . . , m 

(4.4b) 

a'  = j;  a 

“a 

-i 

(4.4c) 

J"'  = r.  r and 

o 

II 

a 

If  the  n standard  deviations  ( Var  (c  ) ^'^^ ) "_  ^ are  given,  then  it 
requires  nXm  divisions  to  compute  all  the  ?^'s,  n multiplications  to 
obtain  a'  and  2n  divisions  to  obtain  r'  and  v'’  • 
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To  facilitate  practical  computations,  it  is  convenient  to  factorize 

C into 

(4.5)  C = Gc’’ 

where  G is  an  m X m )ower  triangular  matrix.  The  factorization  is  pos- 
sible because  C is  symmetric  and  positive  semi-definite.  It  can  be 
achieved  for  example,  by  tlie  scheme  described  in  [7,  p.  491-492].  The 
amount  of  operations  (i.e.  multiplications  and  additions)  required  by  the 
scheme  is  O(m^)  . Since  m is  small,  this  factorization  will,  presumably, 

present  no  computational  difficulty.  With  (4.5),  we  then  compute  the  nXm 
matrix 

(4.6)  V = (p’ p^^)G  . 

The  computation  of  1'  requires  nXm  multiplications  and  additions.  By 
(4.5)  and  (4.6),  (4.3)  can  lie  rewritten  ns 

(4.7)  M = 1 + rr’''  . 

, T 

As  we  shall  see,  there  is  no  need  to  compute  the  product  IT  . The  total 
amount  of  operation;-;  to  Iransforin  the  original  problem  (1.2)  (or  (1.7))  with 
a covariance  matri.-:  givin  by  (4.1)  to  the  equivalent  problem  (4.2)  v/ith  the 

3 

matrix  M given  l.y  (4.7)  is  0(n)  l O(mXn)  l 0 (m‘  ) . 

Recall  that  wc  have  mentioned  some  updating  procedures  in  .Section  2 

which  can  be  used  to  execute  .^tepi  1 in  each  of  the  four  algorithms  in 

2 

Section  3.  These  jirocedures  usually  require  0(n  ) storage  and  are  operat- 
ing on  matrices  with  varying  sizes  which  sometimes  can  be  very  large.  On 
the  other  hand,  by  exploiting  the  special  structure  of  the  matrix  M given 
in  (4.7),  wo  can  achieve  enormous  savings  in  both  the  storage  and  opera- 
tions of  the  algorithms.  More  specifically,  we  can  reduce  total  storage  to 

2 

0(n)  + O(inXn)  + 0(m  ) and  in  .Step  1 we  need  to  operate  on  matrices  of  order 
m only.  These  arc  indeed  savings  because  m is  usually  vci'y  small  com- 
pared with  n . In  what  follow.s,  v;o  illustrate  tiow  those  savings  can  be 
acliieved  in  the  case  witli  upper  liounds  and  omit  the  easier  case  without  the 
upper  bounds. 
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In  the  sequel,  wc  put  subscripts  in  the  identity  matrix  I to  denote 
its  dimensionalities  and  write  I SI  to  moan  the  cardinality  of  a set  S . 
Let  r be  an  n X ni  matrix  and 

nT 


(4.7) 


M = + rr^ 


Let  p,q,r  and  a be  n-vectors.  Let  I,  and  subsets  of  {l,...,n}, 

disjoint  from  each  other  and  such  that  I ti  = {l,...,n}  . Then  wo 
have 

Vi  ^ ''^i'  '^1- 

whose  inverse  can  be  computed  by  the  Sherman-Morrison-Woodbury  formula  ([11, 
p.  79]) 

= I,  ^ I - [I  + (J\  )'^()'  )]“^  (r^  )"^  . 

I J I J 

For  each  6 rR  ^ , let  a be  the  (unique)  solution  of 


(4.8) 


[I  + (I\  )'^n\  )lo*  = (r_  )'^  <>;  . 

In  ^1* 


Now  if  f-  , g.  and  h^  are  the  solutions  of  tlic  systems  of  equations  in 

Jf 

(3.20) , then 

Pj 

(4.9)  f^=r-r,  a ^ and  h_  =Pt  o ^ 

•^1  '^1  -^l.  -^1  >^1  ‘^i- 

To  compute  g.  , it  is  necessary  to  evaluate  M ^ M _ a . Ke  have 

Jl  "^1"^1  '^1"’2  ^2 

= ^Jl/^J2.>%2  ■ 


Thus , 

lllz/:  1*  2 • 2 1 1*  1*  1*  .1*  2.'  I 


= r,  Ii“{i+(r^  )'^(r,  ))~^(r^  )’'(r.,  )](r.  )‘a, 

J,  m m j,  j,  j,  j-  j,. 

1*  1*  1*  1' 
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( 4.13b) 


(4.14a) 


+ ^PjJ  -rj  (o  ^ i - t 2,  ^ 


(4.4b) 


rj  +11  Pj 

■"  PPjJ  - Tj,  (o  ^ 2' 


2 “2  “2  “2 
If  we  define  the  n-voctors  s and  t by 


(4.15a) 


(4.15b) 


s = (q  + Xp)  - r(  0 ^ - r 


t ^ (r  4 pp)  ~ r(a 


then  we  can  iinmcdintely  generate  all  the  vectors  ,t^  ,c'  ,d_,c,  and  d 

Jl  J-i  I 1J-,  k. 


"■^1  '^1  "^1 

In  practice  the  vectors  o ,a  ,a 
the  following  three  systems  of  equations 


and  T are  generated  by  solving 


(4.16) 


‘I  I « T Pi 


= ((!'  )^  a -(r  jr  )V  4 (Tj  Pj  ) • 

Jl.  J2.  J2'  Jj. 

'^1  “2 

Note  that  there  : c no  need  to  generate  0 and  t separately.  It  is 

• T 

their  diffeience  that  is  required.  I'urtliermore , the  matrix  I + (f^  ) (I\  ) 

in  j 2 2 

is  symmetric,  positive  definite  and  of  order  m which  is  independent  of 

1 Jj^l  . Now,  suppose  ttiat  changes  by  one  element,  tlien  it  is  easy  to 

T 

see  that  I + (I'_  ) ( ) isj  affected  by  a i‘ank-onc  modification,  i.c. 

Ill  j 2 *J  2 


(4.17) 


new  matrix  = old  matrix  t yy ' 


wlicre  y is  a suitable  row  in  the  matrix  P . If  wc  solve  (4.16)  by  iipply- 

T 

ing  Cholos);y  factorization  (see  (11]  c.g.)  on  the  matrix  I + (P.,  ) (P^  ), 

■ m cl  , v)  m 

1 • 1 • 

Uicn  we  can  update  the  Cholos);y  factors  very  quickly  by  taking  advantage  of 
the  relation  (4.17).  An  updating  procedure  is  described  in  (8,  Algorit)im  cl ] . 
See  also  [91.  The  number  of  operations  necessary  to  compute  the  modified 
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iactori zntion  using  the  procedure  is  m + 0(m)  multiplications  and 
2 

m + 0(m)  additions  [8].  Furthermore,  the  procedure  is  numerically  stable 

T 

because  for  every  J, , the  matrix  I + (1,  ) (F^  ) is  sufficiently 

1 jn  J , 

positive  definite,  its  smallest  eigenvalue  is  greater  than  or  equal  to  1 . 

We  note  that  if  m = 1,  (corresponding  to  the  single  index  model),  then 
T 

+ (F  ) (F  ) becomes  a scalar  and  the  updating  process  reduces  to 
rn  ^ 3.  • i • 

trivial . 

Based  on  the  above  analysis,  we  may  now  formulate  the  specializations 
of  Algor itlims  6 and  7 to  problem  (4.2)  with  the  upper  bounds.  To  keep  the 
notations  simple,  we  delete  the  primes  in  all  the  vectors  in  (4.2).  More- 
over wc  assume  that  the  vector  q does  not  vanish  in  the  objective  function 

It  t 

of  (4.2),  i.e.  we  consider  the  objective  function  to  be  ^ z Mz  (q  + Or)  z 

T 

Furtheriiiore , we  assume  p at  1 in  order  to  eliminate  the  trivial  case. 
Finally,  we  3'emark  that  the  vectors  g.^  ,h  ,c  ,d  ,c  and  d required  in 
Algoritlim  6 can  be  obtained  easily  from  (4.15)  bv  some  simple  modifications. 
Stn  (4.18)  and  (4.19)  below. 


Algoritlim  8.  Specialization  of  Algorithm  6 with 

Step  0.  Determine  X = max  {-q./p.}  and  lot 

^ l<i<n  " ^ 

Let  l = {l,...,n}\{k),  Jj^=Tk)  and  = <>  ■ 


M given  by  (4.7). 

k be  a maximizing  index. 
Go  to  Step  1 . 


Step  1 . Solve  the  systems  of  linear  equation"? 


(4.18) 


and  compute 


II 


) ] (o 


(4.19) 
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Determine 


(4.20) 


s . 

X = max{max{-  : t.  >0,  irJ),  max{-  r : t.  <0,  j,  cJ,  } , 

new  ^ ^ 

1 3l  A 


a . +s . 

3l 


a , +s  . 

32  32 


mnx{-  — i i : t.  >0,  j cJ,l,  max{-  — : t.  <0,j,cj,)} 

3i  1 J 32  ^ ^ 


Case  (i)  ♦ (t' . Vwt 


— T T T 

A = -n  + .s  -pa  )/p  t 


^ " '^new'^olcl’'  set  A = A 

x‘*(0)  = 0,  (0)  ^ -s,  - A*t^  and  x^  (0)  = a, 

1 J2  J2 

and  terminate.  Ot  lii-rw  i pe , lot  k be  a maximizing  index  in  (4.20).  Set 

A - A and  cio  to  Step  3. 

old  new  ■ ‘ 

Case  (ii)  . If  = 3.  set  A*  = . 

x‘^(0)  0 and  x^  (0)  = a, 

^ ^2  ^2 

and  terminate.  Otlierwiso,  lot  k be  a maximizing  in  (4.20).  Set  = 

A and  qo  to  Step  3. 

new  ^ 

Stop  3 . If  k I I , replace  T and  J by  J\  {k)  and  u {k}  respective- 
ly. Go  to  Step  1 . 

If  krd^/  replace  and  by  J^u  {k}  and  J2^  respectively 

Go  to  Step  1 . 

If  k c and  tj^  < 0 , replace  I and  by  IU{k)  and  \ (k) 

respectively.  Go  to  Step  1. 

If  k f and  replace  and  by  J^\  {k}  and  J2U{k) 

respecively.  Go  to  Step  1 . 

Algorithm  9:  Specialization  of  Alcjoritltm  7 v.’ith  M given  by  (4.7). 

l^e£_0.  Use  Algori'hm  S to  obtain  x^(0)  and  let  I,  and  be  the 
three  index  .sct.s  at  t rmination  of  the ! algorithm.  Let  = 0 , If 


> J 


* 

= <t>  set  s = q + A p anc3  t = r,  yo  to  Step  2.  Otiiorwiso  go  to  Stop  J 
Stop  1.  Solve  tlie  systems  of  equations  in  (4.16).  Compute 


T T 
n = P,  1 , 

Jl  Jj. 

Pf.  _ , T ^ T,  ^1  ^2. 

PG  - p q - n (ii  - X ) 


r, 

T T 1 

PF  = pi,  r - ,10 

P, 

T T 1 

PH  = p'  p - „'o 


X = - n H PG  - p"  a )/I’li,  u = -PF/PH 
‘ 2 ‘^2 


and  the  vectors  s and  t by  (4.16). 

Step  2.  If  tT>0,  t^  “ h , I <0  -not  -=  0,  ((,)  - -s_  and 

I-  I Jj 

x‘  (0)  -=  a_  for  0>0  cind  lenuinato.  Otliorwasc  doVonnino 

J 2 a 2 **  o I d 


R . 3 . 

= min  fmin  {-  t--:  1. . < 0 , i f I } , min  {-  r : t . > 0 , j , c J,  ) , 

t.  a t.  3,  '-^1  1 


3, 


lin  f- 


a . +s  . 

Ij  Ij 


a . Hr. 
^2  •‘’2 


f t <0,  j^t  mini-  t 

^1  ^ 32 


0,  io  r d., } } 


and  let  k 
'a 


be  a minimising  index.  Put  x'  (0)  = 0,  >:  _ (6)  = -s.,  - 

X t/  • O « *J  ^ 


t.^  and 

3 '^l  1 

x“  (0)  -■  a,  for  Or  lO  , ,,0  ] . Set  0 , , = 0 and  no  to  Stop  2. 

old  nc'V  old  new  ■'  ‘ 


Step  3.  Same  ns  Step  3 in  AlqoriUim  7. 

Tlic  systems  of  linear  equations  in  (4.16)  or  (4.18)  should  best  be 
solved  by  Cliolcsky  factor!  xa(  icni  togtd  lier  will)  Algorithm  Cl  in  [8]  to  update 
the  Cholesky  factors.  Pccau'ic  of  the  fact  that  starts  initially  in 
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Algorithm  8 as  n singleton,  there  is  no  need  to  compute  t}io  product 
T 

(r,  ) (I't  ) throiinhout.  Final]  whenever  J,  becomes  a singleton,  it 

•^1.  ^ 

is  more  convenient  to  compute  the  Cholesky  factors  of  the  matrix 

T 

1 + (r,  ) (1-,  ) from  scratch.  Formulas  for  computing  such  factors  arc 

ni  J ^ J 2^ 

given  in  [8]. 

Of  course,  the  above  two  algorithms  are  applicable  to  the  general 
portfolio  analysis  problem  (1.7)  with  the  covariance  matrix  V given  by 

(4.1) ,  provided  that  wc  first  obtain  the  transformed  problem  (4.2)  with  the 
transf  oj  .mod  data  (4.3)  and  (4.4)  liy  the  nrocess  described  earlier. 

The  total  storage  required  in  tliis  anpl  i cat;- on  consi.sts  of  the  n standard 
doviatic;ns  Vfir  (te  ) ^^^ , the  mXm  covariance  matrix  C,  the  mXm  matrix 
G,  (C  and  G each  require  storage  because  of  symmetry),  the  nxm 

matrices  ^ Py ' ’ ■ ‘ n-vcctors  r,a  and  x,  finally,  some 

storage  for  the  Cholesky  factors.  These  add  up  to  0(n)  -t  O(m^) 

4 O(nXni)  . This  latter  amount  is  best  described  as  minimal. 

We  concl  iid(.-  tliis  .■'ection  liy  remarking  that  the  matrix  V given  by 

(4.1)  appeared  also  in  an  application  of  the  parametric  linear  complemontar 
ity  problem  to  structural  engincei-ing.  Jn  fact,  it  was  shown  in  115)  (see 
also  the  references  tliere)  that  the  behavior  of  a rciinforccd  concrete  beam 
can  be  described  in  terms  of  the  iiarnmetr  ie  linear  eoraplcment  a I'i  ty  problem 

{(q  -I  ip>,M):  A r U}  where  M is  precisely  c>f  the  form  (4.1).  Of  course. 
Algorithm  2 is  a]>pl icable  in  this  insLaneo.  Moreover,  the  analysis  devel- 
oped in  this  section  c^in  be  used  to  increa.se-  the  overall  efficiency  of  the 
algorithm.  Pee  [23]  for  a more  detailed  di:;cussion  of  this  subject. 


5.  CONCLUSION 


In  this  paper,  wo  have  developed  several  new  algorithms  for  portfolio 
analysis  and  discussed  their  efficient  implementation.  We  have  also  shown 
how  their  specializations  to  the  index  models  will  result  in  dramatic 
savings  in  storage  and  computations.  In  a forthcoming  paper  (23],  we  shall 
report  our  computational  experience  with  the  algorithms  in  solving  some  port- 
folio problems. 
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